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THE  LAST  14  DAYS  OF  SKYLAB  1:  ORBIT  DETERMINATION  AND  ANALYSIS 

by 

Doreen  M.  C.  Walker 


SUMMARY  '  ^ 

The  orbit  of  Skylab  1  (1973-27A)  has  been  determined  using  some  1400 
NORAD  observations  during  the  14  days  prior  to  decay  on  1979  July  11.  There 
are  14  daily  orbits,  with  standard  deviations  corresponding  to  average 
accuracies  of  40  m  cross  track  and  30  m  radial.  A  15th  orbit,  only  slightly 
less  accurate,  was  determined  from  observations  on  July  11  between  the 
manoeuvre  at  07.45  UT  and  decay  at  16.37  UT. 

The  variations  in  inclination  due  to  atmospheric  rotation  and  16th- 
order  resonance  with  the  geopotential  have  been  successfully  analysed,  to 
give  the  first  values  of  16th-order  geopotential  coefficients  determined 
from  resonance,  and  a  value  of  1.10  ±  0.07  rev/day  for  the  atmospheric 
rotation  rate  at  a  height  of  210-220  km. 


The  daily  changes  in  semi-major  axis  have  been  used  to  determine  13 
daily  values  of  air  density,  at  heights  from  252  km  down  to  179  km.  All 
agree  well  with  the  CIRA  1972  model,  and  indicate  a  smaller  semi-annual 
variation  than  in  the  early  1970s. 

The  variations  of  eccentricity  and  argument  of  perigee  take  unusual 
forms,  but  detailed  analysis  shows  that  the  variations  are  in  full  accord 
with  the  theory  for  an  atmosphere  with  day-to-night  variation  in  density, 
with  the  perigee  progressing  towards  the  point  of  minimum  density., 
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INTRODUCTION 


Skylab  1,  1973-27A,  was  launched  on  14  May  1973  at  17.30  UT  into  a  nearly  circular 

orbit  inclined  at  50°  to  the  equator*.  The  central  structure  of  the  spacecraft  was 

2 

cylindrical  ,  25.6  m  long  and  6.6  m  in  diameter,  while  its  solar  panels  gave  it  an  over¬ 
all  span  of  27.4  m  (see  Fig  1),  and  its  mass  was  74783  kg. 

Three  separate  crews  manned  Skylab  between  23  May  1973  and  8  February  1974,  and 

after  this  date  the  spacecraft  was  left  to  decay  naturally.  The  total  lifetime  for 

3 

Skylab  was  then  predicted  by  RAE  as  6  years  .  This  prediction  proved  to  be  very  good, 
well  within  the  usual  error  for  lifetime  predictions,  10Z  of  the  remaining  lifetime.  An 
account  of  subsequent  prediction  by  RAE  is  given  in  Ref  4.  Skylab  finally  decayed  on 
11  July  1979  at  16.37  UT  over  SW  Australia  in  the  full  glare  of  press  publicity,  after 
2249  days  in  orbit.  * 

After  the  decay,  over  2000  observations  of  the  satellite,  made  during  its  last 
15  days  in  orbit  by  the  assigned  and  contributing  sensors  of  the  North  American  Air 
Defense  Command  (NORAD),  Space  Detection  and  Tracking  System  (SPADATS),  were  made 
available  to  us,  and  in  this  Report  the  observations  have  been  analysed  to  give  daily 
orbits  using  the  RAE  orbit  refinement  program  PROP  in  the  PROP  6  version3. 

2  THE  ORBITS 

2. 1  Up  to  July  9 

In  the  PROP  6  model  the  mean  anomaly  M  is  c-',-«-ed  by  a  polynomial  of  the  form 

M  -  Mq  +  Mjt  +  M2t2  +  M3t3  +  M4t4  +  M5t5  ,  (1) 

where  t  is  the  time  measured  from  epoch,  and  the  number  of  M  coefficients  used 
depends  on  the  severity  and  variability  of  the  drag.  The  orbital  elements  for  the  first 
13  orbits  determined  daily  from  June  27  to  July  9  are  given  in  Table  1,  with  the 
standard  deviations  below  each  value.  The  epoch  for  each  orbit  is  at  00  h  on  the  day 
indicated. 

As  the  individual  orbits  were  determined  from  observations  extending  over  less 
than  24  hours,  so  that  t  <  0.5  ,  it  would  be  expected  that  only  a  small  number  of 
coefficients  in  equation  (1)  would  be  required.  This  proved  to  be  true,  with  seven 
orbits  needing  only  Mg  -  M2  and  the  remaining  six  Mg  -  M^  .  The  orbit  determined  on 
June  27,  using  only  25  observations  poorly  distributed  around  the  orbit,  was  not 
reliable,  the  value  of  inclination  in  particular  being  badly  in  error.  So  in  the  subse¬ 
quent  analysis  this  orbit  (given  in  square  brackets  in  Table  1)  was  ignored,  the  orbit 
for  June  28  being  regarded  as  orbit  1.  The  orbit  for  June  27  is  included  in  Table  1 
to  illustrate  the  bias  errors  that  can  occur  when  the  distribution  of  observations  is 


Values  of  the  orbital  parameters,  with  standard  deviations,  for  orbits 
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inclination  (deg)  e  =  measure  of  fit 

right  ascension  of  ascending  node  (deg)  N  =  number  of  observations  used 

argument  of  perigee  (deg)  D  =  time  covered  by  the  observations  (days) 
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For  the  12  daily  orbits  from  June  28  to  July  9  (orbits  1  to  12),  the  observations  for 
each  orbit  are  much  more  numerous  and  are  fitted  in  a  satisfactory  manner,  with  e  ,  the 
parameter  indicating  the  measure  of  fit,  always  between  0.24  and  0.40.  It  would  appear 
that  the  range  error  allocated  to  the  NORAD  observations  could  be  reduced,  because  the 
range  residual  is  very  small  for  the  great  majority  of  observations.  For  these  12  orbits 
the  average  standard  deviation  in  inclination  is  0.0003°,  equivalent  to  35  m  in  cross¬ 
track  distance,  and  0.000004  in  eccentricity,  equivalent  to  30  m  in  radial  distance. 

The  last  column  in  Table  1  gives  the  mean  height  of  the  satellite,  6372  km  being 
the  appropriate  mean  Earth  radius. 

2.2  For  July  10.0 

During  its  last  2  days  in  orbit  Skylab  was  tracked  by  nearly  every  available  radar 
because  it  was  known  that  quite  large  sections  would  survive  re-entry  and  fall  somewhere 
on  the  Earth's  surface.  The  idea  was  to  obtain  observations  and  keep  updating  the  impact 
prediction,  and  if  it  seemed  that  Skylab’s  last  orbit  would  be  over  a  highly  populated 
area,  attitude  changes  would  be  made.  In  consequence,  many  hundreds  of  observations  were 
available;  about  600  from  midday  on  July  9  to  midday  on  July  10,  and  about  800  between 
midday  on  July  10  and  decay  at  16.37  UT  on  July  11. 

The  600  observations  centred  on  July  10.0  fell  into  two  groups  of  about  300,  one 
before  and  one  after  midnight.  At  first,  as  PROP  will  only  accept  100  observations,  an 
attempt  was  made  to  determine  orbits  from  observations  extending  over  only  about  3  hours 
(2  revolutions).  However,  these  orbits  proved  to  be  unreliable  because  the  observations 
were  concentrated  at  a  small  number  of  stations  with  inadequate  geographical  distribution, 
leading  to  geometrical  bias.  Also  the  orbital  period  and  its  rate  of  change  were  ill- 
defined,  because  the  time  interval  of  the  observations  was  too  short,  and  the  drag  was 
so  high  that  could  not  be  omitted  from  the  model. 

So  the  time  interval  was  extended  to  12  hours  -  before  and  after  midnight.  This 
meant  that  two  thirds  of  the  observations  had  to  be  omitted,  and  this  was  accomplished  by 
reducing  the  number  of  observations  per  transit  at  each  station.  At  some  stations  there 
were  more  than  20  observations  per  transit  so  the  removal  of  two  thirds  was  not  serious. 

The  two  sets  of  elements  obtained  are  given  in  Table  2  as  13A  and  13B  with  the  standard 
deviations  below  each  value:  the  first  orbit,  13A,  is  determined  from  observations  before 
midnight  (July  10.0)  and  the  second,  13B,  from  observations  after  midnight.  A  third 
orbit  was  determined  from  observations  taken  from  both  13A  and  I3B  at  times  within  6  hours 
of  the  epoch,  July  10.0.  This  orbit  is  given  as  orbit  13C  in  Table  2. 

The  agreement  between  the  three  orbits  at  epoch  July  10. 0  is  very  good.  The  values 
of  eccentricity,  e  ,  and  argument  of  perigee,  u  ,  agree  to  within  their  standard 
deviations,  the  three  values  of  inclination,  i  ,  are  within  twice  the  sum  of  their 
standard  deviations,  and  the  values  of  right  ascension  of  the  node,  0  ,  within  three 
times  the  sum  of  their  standard  deviations.  The  values  of  Mj  (and  hence  semi  major 
axis,  a)  and  M2  however  do  not  agree  so  well.  This  is  because  orbits  I3A  and  13B  are 
determined  from  observations  away  from  epoch  and  therefore  experience  different  drag 
conditions.  Since  the  model  (for  orbit  13)  does  not  include  M3  ,  the  values  of  M2  in 
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later  coefficients  in  the  polynomial  for 


orbits  I3A  and  13B  are  those  appropriate  over  the  time  interval  of  the  observations, 
centred  at  July  9.7  and  July  10.3  respectively.  The  value  of  M2  for  orbit  13C  should 
be  that  appropriate  to  July  10.0,  being  derived  from  observations  over  the  epoch,  and 
the  value  of  on  orbit  13C  should  be  taken  as  correct,  since  the  values  on  orbits 

13A  and  13B  are  extropolated  using  values  of  M2  derived  from  the  fitting  of  observa¬ 
tions  away  from  epoch. 

2.3  For  July  1 1 .0 

On  the  morning  of  July  11  at  07.45  UT  Skylab  was  commanded  to  manoeuvre  to  a 
tumble  attitude,  thus  reducing  the  drag  and  extending  the  lifetime.  This  action  was 
taken  to  shift  the  probable  impact  footprint  away  from  the  highly  populated  east  coast 
of  the  United  States  and  Canada  to  the  Indian  Ocean0.  Therefore  the  800  observations 
available  between  midday  on  July  10  and  decay  were  divided  into  three  groups:  about  400 
between  12.00  h  and  midnight  on  July  10;  about  90  between  00  h  and  the  manoeuvre  at 
07.45  UT  on  July  11;  and  finally  some  300  observations  between  07.45  UT  and  decay. 

For  the  orbits  derived  from  the  first  two  sets  of  observations,  14A  and  14B,  at 
epoch  July  11.0,  it  was  necessary  to  alter  the  orbital  model,  because  the  eccentricity 
increased  slightly  between  July  9.0  and  July  11.0  (due  to  the  day-to-night  variation  in 
air  density),  whereas  the  PROP  model  indicates  a  decrease  in  accordance  with  the  theory 
for  decay  in  a  spherically  symmetrical  atmosphere.  At  July  11.0  the  rate  of  change  of 
e  per  day  given  by  the  model  was  strongly  negative  (-0.00011  for  14A  and  -0.00016  for 
14B) ,  and  an  appropriate  positive  e  was  added  to  the  model  to  give  consistent  results. 
The  method  chosen  was  to  adjust  the  added  5  so  as  to  give  identical  values  of  e  on 
the  two  orbits:  the  added  4  was  0.00010  on  both,  giving  4  on  the  (modified)  model  as 
-0.00002  for  14A  and  -0.00005  for  14B.  (These  values  are  different  because  the  values 
of  M2  are  different.)  This  procedure  adopted  for  14A  and  14B  was  not  necessary  on 

orbits  13A  and  13B,  presumably  because  the  rate  of  decrease  of  e  in  the  PROP  model  was 

much  smaller  (being  proportional  to  M2). 

The  sets  of  elements  for  orbits  14A  and!4B  are  given  in  Table  2.  A  third  orbit  was 
determined  from  observations,  taktn  from  14A  and  14B,  over  the  epoch,  July  11.0.  For 
this  orbit,  14C,  the  observations  ranged  from  17.56  UT  on  July  10  to  04.23  on  July  11. 

The  orbital  elements  for  orbit  14C  are  given  in  Table  2.  Comparing  the  three  sets  of 

elements  for  July  11.0,  orbits  14A,  B  and  C,  we  see  that  the  values  of  Mj ,  M2  and  M^  , 

and  hence  a  ,  are  much  more  accurately  determined  for  orbit  I4C,  and  this  is  probably 
due  to  the  wider  spread  of  the  observations  in  time  and  perhaps  in  latitude,  and  the 
fact  that  they  cover  the  epoch.  Comparing  the  values  of  e,  i  and  0  on  orbit  I4C 
with  those  on  orbits  I4A  and  14B  it  can  be  seen  that  they  agree  to  within  2.3  times  the 
sum  of  their  standard  deviations.  The  values  of  M2  on  orbits  14A  and  14B  differ  and, 
as  before,  may  be  taken  as  applying  at  a  time  midway  through  the  observations  rather  than 
at  July  11.0. 

In  the  analysis  to  follow,  the  orbits  I3C  and  I4C  have  been  used  for  July  10.0  and 
11.0  respectively. 
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2.4  For  the  last  eight  hours  (08.08  to  16.02  UT)  on  July  11 

The  observations  available  after  the  tumble  action  at  07.45  UT  on  July  II  have  been 
used  to  derive  one  further  orbit.  As  the  epoch  at  which  the  orbit  is  determined  by  the 
PROP  program  is  always  at  00  h,  the  observations  were  being  fitted  away  from  epoch.  So, 
besides  giving  the  orbit  determined  at  00  h  as  14D  in  Table  2,  the  elements  have  been 
converted  to  midday  to  give  an  orbit  15  for  July  11.5.  Orbit  15,  with  an  epoch  in  the 
midst  of  the  observations,  is  obviously  preferable  for  use  in  any  further  analysis. 

There  was  a  difficulty  with  the  convergence  process  in  the  PROP  program  for  this 
last  orbit  (14D),  and  results  could  only  be  obtained  by  fixing  the  value  of  e  and 
allowing  the  other  parameters  to  change.  After  this  another  PROP  fitting  was  run  keeping 
all  the  values  determined  on  the  previous  run  fixed,  but  leaving  e  free  to  be  deter¬ 
mined.  However,  the  change  in  e  was  negligible,  only  1  *  10  0  .  There  was  still 
the  problem  of  estimating  the  value  of  the  standard  deviation  for  e  .  This  was  done  by 
comparing  the  standard  deviations  printed  out  on  runs  which  had  not  converged  but  had 
reached  the  same  value  of  e  as  the  run  with  e  fixed.  On  this  basis  the  standard 
deviation  for  e  was  estimated  as  10  *  10  ®  . 

Although  orbit  14D  is  at  an  epoch  well  outside  the  time  span  covered  by  the 
observations,  the  values  of  e,  i,  0  and  w  are  surprisingly  good.  The  values  of 
Mq,  M,  (and  a  )  and  on  orbit  14D  look  peculiar,  but  this  is  merely  because  they 

are  coefficients  of  a  polynomial  with  an  inappropriate  zero  point,  and,  when  converted 
to  epoch  July  11.5,  they  give  reasonable  values.  These  values  in  orbit  14D  have  been 
enclosed  in  square  brackets  to  emphasize  their  inappropriate  character. 

The  accuracy  of  orbit  15  cannot  be  formally  assessed  because  of  the  convergence 
problem  on  14D,  but  since  a  good  value  of  e  was  achieved  (only  40%  higher  on  I4D  than 
on  14C),  it  would  be  surprising  if  the  main  orbital  parameters  (.ie  excluding  Mq,  M(  and 
M£>  had  errors  more  than  twice  as  large  as  14C.  However,  the  errors  are  likely  to  be 
larger  than  on  14C.  Thus  the  standard  deviations  on  140  seem  suitable  for  use  on 
orbit  15. 

3  ANALYSIS  OF  THE  INCLINATION 
3. 1  Treatment  of  the  data 

The  values  of  inclination  on  the  15  orbits  from  Tables  1  and  2,  with  their 
standard  deviations,  are  plotted  in  Fig  2.  The  values  for  orbits  13  and  14  are  taken 
from  the  orbits  from  observations  over  the  epoch  -  namely  13C  and  I4C.  The  general 
decrease  due  to  the  effect  of  atmospheric  rotation  is  visible,  but  all  other  perturbations 
must  be  removed  before  an  analysis  can  be  attempted  and  a  value  of  the  upper-atmosphere 
rotation  rate,  A  ,  determined. 

Fig  3  gives  the  values  of  inclination  plotted  against  date  after  removal  of  the 

zonal  harmonic,  J,  .  and  lunisolar  perturbations.  The  zonal  harmonic  and  lunisolar 
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perturbations  were  removed  by  using  the  PROD  computer  program  with  1-day  integration 

steps,  and  the  J^  ^  tesseral  harmonic  perturbation  by  using  the  value  recorded  with 

the  PROP  printout. 
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During  Che  14  days  of  the  analysis  the  perturbation  m  inclination  due  to  earth  and 
ocean  tides®  should  not  build  up  to  more  than  0.0001°  (10  metres)  and  can  therefore  be 
absorbed  in  any  standard  deviations  of  0.0003°  or  greater;  however  the  three  standard 
deviations  of  0.0002°  are  increased  to  0.0003°,  to  allow  for  neglect  of  earth  and  ocean 
tide  effects,  and  this  increase  is  shown  in  Fig  3.  Similarly  the  estimated  change  in 
inclination  in  14  days  due  to  solar  radiation  pressure  is  less  than  0.0001  ,  and  can 
also  be  absorbed  in  the  standard  deviation. 

There  remains  one  source  of  perturbations  which  cannot  be  adequately  modelled. 
Because  of  the  peculiar  shape  of  Skylab  1,  Fig  1,  there  is  a  possibility  of  aerodynamic 
forces  perpendicular  to  the  orbital  plane,  which  may  perturb  the  inclination  enough  to 
degrade  the  accuracy  of  the  determination  of  wind  speed.  This  is  discussed  in 
section  3.5. 


3.2  Analysis  of  the  variation  of  inclination  with  time 

The  theoretical  variation  of  inclination  due  to  atmospheric  rotation  and  meridional 
winds  can  be  calculated  for  a  series  of  values  of  A  ,  and  of  p  ,  the  south-to-north 
atmospheric  rotation  rate,  using  a  computer  program  (ROTATM)  based  on  equations  (32)  and 
(33)  of  Ref  10. 


The  variation  was  computed  for  values  of  A  from  1.0  to  1.35  rev/day  at  intervals 
of  0.05  with  u  ■  0  and  0.1  rev/day.  The  addition  of  the  u  -  0.1  term  had  no  effect. 


This  was  not  surprising  because  a  constant  south-to-north  wind  has  no  effect  on  the 
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inclination  of  a  circular  orbit  .  Even  with  a  non-circular  orbit,  the  variation  in 
inclination  due  to  meridional  rotation  is  proportional  to  cos  u  ,  and  here  <u  is  near 


90°  until  after  orbit  10  on  July  7;  during  the  last  4  days  w  changes  from  53°  to  332° 


and  the  local  time  is  04  ±  3  h,  when  the  meridional  wind  is  normally  towards  the  equator, 
and  as  perigee  moves  from  north  to  south  of  the  equator  the  effect  of  the  wind  would  tend 
to  cancel  out.  So  there  is  no  reason  to  suppose  that  meridional  winds  are  of  any 
importance  here. 


The  values  of  inclination  were  therefore  fitted  by  choosing  the  best  vdlue  of  A  , 
assuming  p  *  0  .  The  best  fit  between  theory  and  the  observational  values  of  inclina¬ 
tion  was  difficult  to  decide,  but  appears  to  be  with  A  between  1.20  and  1.30.  The 
theoretical  curve  for  A  ■  1.25  has  therefore  been  drawn  through  the  points  in  Fig  3. 

In  order  to  assess  the  likely  errors  in  the  values  of  A  ,  a  realistic  estimate  of  the 


likely  errors  in  i  at  the  beginning  (i„)  and  end  (i„)  of  the  curve.  a„  and  o„  ,  was 

™  '  1  D  b  DC* 


n  i  /.  * 

made,  and  then  KJ  +  o^  /  (ig  -  ig)  was  taken  as  the  standard  deviation  in  A  .  This 
was  found  to  be  0.11. 


The  value  of  A  ■  1.25  i  0.U  obtainec  from  fitting  the  values  of  i  in  Fig  3,  ca 
be  regarded  as  averaged  in  local  be'  -e  the  orbit  is  so  nearly  circular.  The 

height  decreases  from  253  km  to  I4<'  ..  t  during  the  13j  days,  and  the  average  height  at 
which  A  applies  is  difficult  to  determine.  In  the  Appendix  an  attempt  is  made  to 
derive  as  logical  as  possible  a  value  for  the  mean  height,  based  on  averaging  the  air 
density.  This  method  gives  a  mean  height  of  211  km,  but  the  omission  of  orbit  15  alters 
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the  mean  height  to  218  km,  while  still  yielding  the  same  value  of  A  .  So  the  mean 
height  should  perhaps  be  regarded  as  a  height  band  of  210  to  220  km,  rather  than  an 
exact  value. 

The  value  of  A  is  larger  than  might  be  expected  from  Fig  13  of  Ref  12,  where  a 
value  of  1.1  for  A  is  indicated  for  a  satellite  at  210  to  220  km  height  in  'average* 
conditions. 

3.3  Analysis  of  inclination  in  terms  of  orbital  period 

For  satellites  in  near-circular  orbits,  with  eccentricity  less  than  0.005,  the 
change  Ai  in  inclination  produced  by  a  change  AT^  in  orbital  period  can  be  expressed 
much  more  simply  than  in  the  lengthy  theoretical  equation  of  Ref  10.  The  simplified  form 

4  ■  tM'  * lo  *  $)]  ■ 

where  z  =  ae/H  is  of  order  0.1.  The  parameter  c  allows  for  the  effect  of 
atmospheric  oblateness  and 

2 

e'a(l  -  e)  sin  i 
2H 

where  e'  is  the  ellipticity  of  the  atmosphere  (taken  as  0.00335)  and  H  is  the  density 
scale  height.  The  factor  F*  in  equation  (2)  is  given  by  ■  1  -  [a(l  -  e)w  cos  i}/Vp 

where  is  the  satellite's  velocity  at  perigee  and  w  is  the  angular  velocity  of  the 

atmosphere  near  perigee.  So,  if  A  is  constant,  the  values  of  inclination  for  a 
circular-orbit  satellite,  when  plotted  against  orbital  period,  should  lie  on  a  straight 
line  whose  slope  gives  the  value  of  A  . 

Here  the  eccentricity  of  Skylab  is  below  0.001 ,  so  this  method  should  yield  a  good 
determination  of  A  .  The  values  of  inclination,  with  perturbations  removed,  are 
plotted  against  anomalistic  period  in  Fig  4  and  a  least-squares  fitted  straight  line  to 
the  values  of  inclination  has  a  slope  of  0.0074  ±  0.0005  deg/min.  For  this  satellite, 
with  c  ■  0.186  and  F'  *  0.908  ,  equation  (2)  gives  A  ■  1.28  ±  0.09  .  This  is 
consistent  with  the  value  found  in  section  3.2,  namely  A  -  1.25  ±0.11  .  But  the  fitting 
is  not  as  good  as  might  be  hoped,  tha  measure  of  fit  e  being  2.5. 

3.4  16th-order  resonance 

As  the  values  of  A  obtained  in  sections  3.2  and  3.3  are  larger  than  might  have 
been  expected,  and  the  fitting  fell  short  of  expectations,  it  was  natural  to  look  for 
any  other  perturbation  that  might  affect  the  value*  of  inclination.  With  this  in  mind, 
the  l6th-order  resonance  was  investigated.  Previously  it  had  been  assumed  that  the 
satellite  would  pass  through  16th-order  resonance  very  quickly  and  that  no  change  in 
inclination  would  be  discernible.  The  date  of  resonance  was  July  8.9,  three  days  before 
decay.  In  the  past  we  have  normally  limited  the  range  of  $  ,  the  rate  of  change  of  the 
resonant  angle,  to  ±20  deg/day  when  the  orbits  available  over  the  resonant  period  were 
7  days  apart.  Here  the  analysis  can  be  extended  to  a  much  larger  value  of  5  ,  as  the 
orbits  are  daily  and  the  range  of  *  over  which  the  resonance  is  analysed  will  be 
- M —  


The  theory  for  the  rate  of  change  of  inclination  at  resonance  is  given  in 
section  5.1  of  Ref  10  where  the  parameters  used  are  defined.  The  first  term,  ie  the 
(y,q)  ■  0  term,  in  the  equation  for  rate  of  change  of  inclination  near  16th-order 
resonance  is 


,,  n  /R\17  -  L°. 1 

-j-  =  — : - :(— )  (16  -  cos  i)F  q)  S  -  sin 

dt  sin  i\a/  17,16,8 \  16 


_0>1 


$  +  Cjg  cos  $ 


(3) 


where  $  °  u  +  M  +  16 (ft  -  v)  is  the  resonance  angle, 


being  the  sidereal  angle. 


The  values  of  inclination,  cleared  of  zonal  harmonic,  J 


2,2 


and  lunisolar 


perturbations,  were  fitted  with  equation  (3),  in  integrated  form,  using  the  THROE 
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computer  program  .  The  first  fitting  was  made  with 


A  »  1.20  and  e  ,  the  measure-of- 
fit  parameter,  was  2.36.  Three  of  the  worst-fitting  values,  with  weighted  residuals  over 
3.0,  were  then  degraded  by  a  factor  of  2  and  the  resulting  lumped  coefficients  were: 


9-°.‘ 
10  C, 6 


149  ±  44 


9- 

10  S 


0,1 

16 


-  34  ±  37 


with  e  »  1 .65  . 

The  values  of  inclination  could  now  also  be  cleared  of  the  resonance  perturbation 
and  a  straight  line  was  fitted  to  the  values  by  least  squares,  as  in  section  3.3.  This 
straight-line  fit  gave  a  value  of  A  •  1.12  ±  0.05  .  Ac  this  was  considerably  different 
from  the  value  of  A  used  in  the  THROE  fitting,  the  program  was  re-run  with  A  »  1.10  . 
This  yielded  the  following  values  of  lumped  coefficients: 

q_0,l  o_0»l 

lO^C,,  -  147  +  42  10  S,,  -  4  ±  35 
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with  e  *  1 .56  after  one  further  value  was  degraded  by  a  factor  of  2  to  bring  all  the 
weighted  residuals  below  2.5.  The  fitted  curve  is  plotted  in  Fig  5. 

The  values  of  inclination  were  again  cleared  of  the  resonance  perturbation  using 
this  revised  fitting  of  equation  (3)  with  A  *  1.10  and  the  resulting  values  fitted 
with  a  least-squares  straight  line.  The  two  values  of  inclination  for  July  10  and 
July  II  had  their  standard  deviations  increased  for  this  fitting  by  a  factor  of  2.  The 
resultant  value  of  A  was  1.10  ±  0.07,  and  the  straight-line  fit  to  the  values  of 
inclination  is  shown  in  Fig  6. 

So  the  value  of  A  obtained  from  the  straight-line  fit  to  the  values  of 
inclination  is  now  consistent  with  the  value  used  in  the  THROE  fitting.  The  higher 
values  of  A  obtained  in  sections  3.2  and  3.3  occur  because  there  is  a  decrease  in 
inclination  due  to  resonance  of  approximately  0.002°,  as  well  as  the  decrease  in 
inclination  due  to  atmospheric  rotation. 

3.5  Discussion 


The  value  of  A  obtained  from  the  straight-line  fitting  to  the  values  of 
inclination,  after  allowing  for  the  resonance,  is  obviously  preferable.  However,  the 
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values  of  inclination  in  Fig  6  do  not  fit  the  straight  line  quite  as  well  as  might  be 
expected,  and  the  parameter  e  indicating  the  measure  of  fit  had  a  value  of  1.4  despite 
the  relaxation  of  a  number  of  values.  The  high  value  of  t  probably  means  that  the 
standard  deviations  obtained  for  the  inclination  from  the  PROP  runs  are  too  low  -  they 
are  equivalent  to  about  40  m.  Previous  values  of  inclination  for  other  satellites 
evaluated  from  similar  observations  have  yielded  standard  deviations  about  twice  as 
large1®’14’13.  The  low  standard  deviations  probably  arise  as  a  result  of  bias  and  over¬ 
fitting  to  long  runs  of  observations  from  the  same  station,  and  60  to  70  m  is  probably  a 
more  realistic  figure  than  40  m.  An  extreme  example  of  such  bias  is  the  orbit  for 
June  27  in  Table  1 ,  which  was  subsequently  ignored  because  it  was  based  on  only  25 
observations  from  only  five  stations.  In  this  orbit  the  value  of  inclination  appears 
to  have  a  bias  error  of  about  0.0067°,  compared  with  a  standard  deviation  of  0.0038°. 


Because  of  the  unusual  shape  of  Skylab,  the  possibility  arises  of  aerodynamic 
forces  perpendicular  to  the  orbital  plane;  but  these  would  probably  be  averaged  over 
each  revolution  and  would  not  then  affect  these  orbits  determined  over  several  revolu¬ 
tions.  As  the  fitting  of  the  line  in  Fig  6  is  not  perfect,  it  is  possible  that  the 
lateral  aerodynamic  forces  had  an  appreciable  effect  causing  a  decrease  in  inclination 
at  a  perioa  of  about  88.8  min  and  an  increase  at  a  period  near  88.3  min.  However  this 
seems  rather  unlikely,  since  the  spacecraft  attitude  was  maintained  in  a  fixed  mode0 
during  this  time.  The  alternative  explanation  already  given,  of  the  standard  deviations 
being  too  low,  is  therefore  to  be  preferred. 


The  fitting  of  the  variations  in  inclination  at  16th-order  resonance  is  the  first 
successful  analysis  of  its  kind.  Previously  all  attempts  to  analyse  16th-order 
resonance  had  been  thwarted  by  the  high  drag  and  the  satellite's  rapid  passage  through 
resonance.  The  analysis  has  succeeded  with  Skylab  1  because  (a)  the  orbit  is  circular, 
thus  minimizing  drag  for  given  orbital  period,  (b)  the  mass /area  is  large,  and 
(c)  numerous  accurate  observations  were  available,  thus  enabling  accurate  orbits  to  be 
determined. 


It  is  of  interest  to  compare  the  values  of  the  lumped  I6th-order  coefficients 

obtained  from  Skylab  1  with  those  given  by  comprehensive  geopotential  models.  The  best 

of  these  is  the  Goddard  Earth  Model,  GEM  10B  '  3  ,  which  extends  to  order  and  degree  36; 

the  terms  of  degree  l  >  36  can  be  taken  from  GEM  10C,  which  goes  to  degree  and  order 

180.  The  lumped  value  C®!1  ,  for  Skylab  1  at  50.0°  inclination,  can  be  expressed  in 

I  o 

terms  of  the  individual  values,  C.  ,  by  the  following  equation 

jL  m  I  O 


.0,1 

C16 


C1 7, 16  *  3*91C19,16  +  6*70c2 1,16  "  4*73C23,16  “  1 *46C25 , 1 6  +  4* 1 SC27 , 16 


-  0.12C29>16  -  3.26C31fl6  ♦  0.44C33fl6  ♦  2.62^^  -  0.340^,,. 


-  2.I7C39  +  terms  of  degree  >  40  , 


and  similarly  for  S  ,  on  replacing  C  by  S  throughout.  If  the  standard  deviations  of 
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the  individual  coefficients  from  GEM  10B  (and  C)  are  taken  as  3  x  10  ,  the  lumped 

values  are  as  given  below 


from  Skylab 

from  GEM  10B/C 

q.0,1 

10  C16 

147  ±  42 

117  ±  34 

q.0,1 

10S>6 

4  ±  35 

35  ±  34  . 

The  agreement  is  very  satisfactory  and  suggests  the  values  from  both  sources  are 
reliable  to  within  their  quoted  standard  deviations. 

4  .  DETERMINATION  OF  AIR  DENSITY  FROM  DECAY  RATE 

4.1  Theory 

Daily  decay  rates  can  be  obtained  from  the  orbits  in  Tables  1  and  2,  and  the  density 
of  the  atmosphere  at  a  mean  height  can  be  calculated  at  daily  intervals.  The  density  p 
at  a  radial  distance  a  ,  the  'mean'  distance  at  which  the  density  values  apply,  can  be 
obtained  by  substituting  T  *  3n(a/y)^A  in  equation  (7.10)  of  Ref  18  to  give 


Aa/At 

p  «  -  i  1  2,  * 

(ya)*6{l  +  Jc  } 

where  Aa  is  the  change  in  a  in  one  day.  Since  ■  631  .35  s  1  and  At  ■  86400  s 
here. 
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10ap 


Aa 


(4) 


54.556(1) *{l  +  U2} 


where  Aa  and  a  are  in  kilometres,  the  area/mass  parameter  6  is  in  m  /kg  and  p 
3 

is  in  kg/m  .  (The  values  of  Aa  are  given  in  Fig  17.) 

Equation  (4)  includes  the  effect  of  atmospheric  oblateness  but  does  not  allow  for 
the  day-to-night  variation  in  air  density.  If  f  is  the  ratio  of  maximum  daytime 
density  to  minimum  night-time  density  and  F  -  (f  -  1 ) / (f  +  1),  equation  (39)  of  Ref  19 
includes  the  effect  of  day-to-night  variation  in  air  density  and,  with  r^  *  a  ,  gives 


|f  -  -  („«>Spsi0 


['  *  !F(  -  d 


z(P  sin  u  +  Q  cos  w)  +  0(e 


,f2)J  , 


(5) 


IQ  is  the  Bessel  function  of  the  first  kind  and  imaginary  argument  of  degree  n 


where 
and  argument  z 


and  P  and  Q  are  direction  cosines.  Here  F  <  0.22  and  z  <  0.136 


so  that  Iq  may  be  taken  as  1.0  with  error  less  than  0.5%,  and  the  JF  term  in 
equation  (5)  is  less  than  0.015.  Thus  equation  (5)  reduces  to  the  form  appropriate  for 
F  ■  0  and  e  ■  0  ,  with  error  <1.5%,  and  equation  (4)  can  be  used,  as  it  is  the  same  as 
the  reduced  form  of  equation  (5),  except  that  the  effect  of  atmospheric  oblateness, 
expressed  through  the  parameter  c  ,  is  added. 
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The  'mean*  distance  a  at  which  the  values  of  density  apply  is  given  by 


where  aQ  and  are  the  values  of  a  from  Tables  1  and  2  one  day  apart  appropriate 

to  each  daily  value  of  density,  and  6  ■  1 /H  .  If  B(a q  “  a|)  is  small,  equation  (6) 
reduces  to 


For  the  derivation  of  equations  (6)  and  (7),  see  the  Appendix. 

4.2  Results 

Values  of  density  were  determined  at  daily  intervals  using  equation  (4)  with 

appropriate  daily  values  of  3  .  The  value  used  for  5  ,  the  area/mass  parameter,  was 

~3  2 

taken  from  Ref  6  as  7.692  x  10  m  /kg  .  The  values  of  density  p  are  plotted  as 

circles  at  the  top  of  Fig  7,  the  values  of  the  solar  activity  index  SJ0  ^  and  the  daily 

planetary  geomagnetic  index  Ap  being  indicated  below.  The  exospheric  temperature  T^ 

appropriate  for  each  day,  according  to  the  COSPAR  International  Reference  Atmosphere 
20 

1972  ( CIRA  1972),  has  been  calculated  taking  account  of  the  values  of  Sjq  7  and  AP 

for  that  day,  and  the  values  of  density  have  been  standardized  to  a  value  of  exospheric 
temperature  of  1110  K,  by  multiplying  the  density  values  p  by  a  factor  G  ,  given  by 
G  ■  p(l 1 10) /f>(Tm)  ,  where  the  values  of  p(1110)  and  pCT^)  are  obtained  from  the  CIRA 
tables  at  the  appropriate  height.  The  values  of  G  are  plotted  at  the  bottom  of  Fig  7 
and  it  can  be  seen  that  the  adjustment  for  solar  and  geomagnetic  activity  is  as  much  as 
15%  on  some  days.  The  decrease  in  G  between  June  29.5  and  July  4.5  reflects  the 
increase  in  Sjq  j  at  this  time;  the  increase  in  G  between  July  7.5  and  July  10.5  is 
partly  due  to  a  decrease  in  S  n  ,  and  Ap  ,  but  also  partly  because  the  satellite  has 
descended  to  a  height  where  the  effects  of  solar  activity  are  smaller. 

The  adjusted  values  of  density  p*  (•  Gp)  are  plotted  against  height  y  in  Fig  8. 
The  height  is  given  by  y  ■  I  -  R  ,  where  R  is  the  mean  Earth  radius, 

R  ■  R(1  -  J c *  sin  i)  ,  calculated  here  to  be  6372  km.  The  values  of  a  are  obtained 
from  either  equation  (6)  or  (7)  at  the  time  of  each  daily  evaluation  of  p*  . 

The  curves  in  Fig  8  are  plotted  from  CIRA  1972  for  three  values  of  exospheric 
temperature,  900  K,  1000  K  and  1100  K.  The  values  of  p*  ,  which  were  standardized  to  a 
value  of  exospheric  temperature  of  1110  K,  lie  in  a  band  centred  on  T^  -  1025  K  .  This 
reduction  in  Ta  is  to  be  expected  because  all  these  values  of  p*  apply  at  a  time  near 
the  July  semi-annual  minimum  in  density.  In  terms  of  density,  the  observational  values 
are  below  the  CIRA  1110  K  standard  by  an  average  of  15%  for  the  first  eight  values, 
evaluated  between  226  and  25 1  km,  and  an  average  of  9%  for  the  remaining  five  values,  at 
heights  of  179  to  221  km.  The  semi-annual  variation  given  in  CIRA  1972  indicates  a 
decrease  of  14%  at  a  height  of  250  km  on  June  29  and  a  decrease  of  13%  at  a  height  of 
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200  km  on  July  9.  The  values  from  che  analysis  thus  agree  well  with  those  of  CIRA  1972, 

within  1%  for  June  28  to  July  6  and  within  4 7.  for  July  7  to  11.  This  indicates  that  the 

semi-annual  variation  in  July  1979  conforms  to  that  specified  by  the  CIRA  model.  This 

21 

result  is  different  from  that  obtained  for  the  early  1970s  from  analysis  of  1971-106A, 

which  indicated  a  20 7.  reduction  in  density  for  early  July  due  to  the  semi-annual  varia- 

21 

tion.  However,  the  semi-annual  variation  changes  from  year  to  year  ,  and  it  appears 
that  the  variation  in  July  1979  conformed  to  CIRA  1972,  which  is  based  on  results  in  the 
1960s. 

22 

In  a  recent  paper  Eisner  and  Yionoulis  find  that  the  semi-annual  variation  in 

1975  to  1978,  at  heights  of  900  to  1200  km,  has  a  smaller  amplitude  than  that  in  Jacchia's 
23 

1977  model  ,  which  is  very  similar  to  that  of  CIRA  1972.  Their  results,  together  with 
those  obtained  here,  suggest  that  the  amplitude  of  the  semi-annual  variation  in  the  late 
1970s  was  smaller  than  in  the  early  1970s. 

5  ANALYSIS  OF  ECCENTRICITY  AND  ARGUMENT  OF  PERIGEE 

5.1  Variation  of  eccentricity  with  time 

The  values  of  eccentricity  e  for  the  15  orbits  in  Tables  I  and  2  are  plotted  in 
Fig  9.  As  with  the  inclination,  the  values  for  July  10.0  and  11.0  are  from  orbits  13C 
and  14C. 

The  decrease  of  eccentricity  due  to  drag  in  an  atmosphere  with  day-to-night 
variation  in  density,  with  a  day-time  maximum  density  'bulge'  at  14  h  local  time,  is 
given  by  equation  (32)  of  Ref  24  as 

Ae  -  -  k|Ij  +  JFIq  cos  <J>p  +  0(e,  iz2)|  ,  (8) 

where  K  is  a  positive  constant  for  each  epoch  and  4  is  the  bulge-perigee  angle. 

P 

Here  z  is  small,  of  order  0. 1 ,  so  1^  3  1.0  and  I(  3  Jz;  therefore  equation  (8)  may 
be  written  as 

Ae  »  -  iK^z  +  F  cos  4^  +  0(e,  iz2)j  .  (9) 

Between  June  28  and  July  5,  the  value  of  F  is  near  0.2  and  is  between  24°  and 

36°  (see  Fig  12),  so  that  0.81  <  cos  4  <  0.91  ;  hence  F  cos  0  is  slightly  greater 

P  P 

than  z  ,  and  e  should  decrease  much  more  rapidly  as  a  result  of  the  day-to-night 
variation  in  density.  The  broken  curve  in  Fig  9  shows  the  theoretical  decrease  of  e  in 
an  atmosphere  without  day-to-night  variation,  and  the  actual  decrease  is  much  steeper,  as 
predicted  by  the  theory,  between  June  28  and  July  5. 

Between  July  6  and  July  II,  4^  increases  steadily  from  42°  to  178°  (see  Fig  12), 
and  cos  4^  becomes  negative  after  July  8,  thus  making  z  +  F  cos  4^  negative  for 
July  9  to  II.  Therefore  an  increase  in  e  is  to  be  expected  at  this  time,  as  is  seen 
in  Fig  9. 

The  day-to-night  variation  in  density  becomes  very  small  below  190  km,  with 
F  <  0.1  ,  and  so  e  would  be  expect  d  to  decrease  in  the  last  12  hours  of  the  life  in 
accordance  with  the  spherical-atmosp  ».re  theory,  and  the  last  orbit  shows  that  this  is 
beginning  to  occur. 
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5.2 

When  che  eccentricity  is  so  small  (<0.001 )  and  the  drag  so  large,  the  argument  of 
perigee,  u  ,  no  longer  undergoes  its  regular  variation  with  time,  controlled  by  the 
gravitational  effects  of  the  Earth's  oblateness;  instead  the  gravitational  effects  become 
subordinate  to  the  effects  of  the  day-to-night  variation  in  air  density. 

In  the  absence  of  air  drag  Che  variation  of  e  and  u  would  be  a  circular  path  in 
the  (e  cos  <s,  e  sin  u)  plane,  the  centre  of  the  circle  being  on  the  e  sin  u  axis  at  a 

_3 

distance  which,  for  i  -  50.0°  and  a  *  6625  km  ,  is  0.84  k  i0  .In  Fig  10  the 
set  of  points  plotted  as  crosses  shows  the  variation  in  e  that  would  occur  during  the 
14  days  of  the  orbit  determination  in  the  absence  of  drag,  as  given  by  PROD2:  as 
expected,  che  form  is  a  circle  centred  near  the  point  (0,  0.84),  with  small  departures 
from  the  circle  caused  by  lunisolar  perturbations.  The  actual  values  of  e  and  u  , 
also  shown  in  Fig  10,  follow  a  very  different  course,  with  very  large  changes  in  w 
near  the  end  of  the  life,  as  shown  in  Fig  1).  To  discover  the  reason  for  this  unusual 

variation,  it  is  necessary  to  consider  the  theory  for  decaying  near-circular  orbits  in 

.  .  .  19 

an  atmosphere  with  day-to-night  variation  m  density  . 


Variation  of  eccentricity  with  argument  of  perigee 


In  the  theory  it  is  assumed  that  the  maximum  daytime  density  is  at  14  h  local  time 

and  the  minimum  night-time  density  at  02  h  local  time,  both  on  the  equator.  (This  is 

an  approximation  made  in  the  theory;  the  maximum  would  actually  be  north  of  the  equator 

1 9 

in  June  and  July.)  The  theory  is  developed  in  terms  of  z  =  ae/H  (with  constant  H) , 

and  shows  that,  as  decay  approaches,  z  tends  to  approach  the  value  zf  -  j F  cos  $p| 

and  w  tends  to  approach  the  value  Wj,  ■  tan  1 (-  cos  i  tan  M' )  ,  where 

M'  ■  ft  -  L  -  30°  and  L  is  the  solar  longitude.  Also  the  value  of  the  bulge- 

perigee  angle  P  tends  towards  the  value  $  _  ,  where 
P  pt 


-  -  |l  -  sin2i  sin2M'| 


The  values  of 


for  each  day  of  the  orbit  determination  are  shown  in  Fig  12, 


together  with  the  actual  values  of  .  Between  June  28  and  July  4,  the  impression  is 

that  $  is  making  no  attempt  to  approach  <p  _  ,  and  this  is  not  surprising  since  the 
P  p* 

drag  is  not  severe  enough  at  this  stage  to  cause  the  perigee  to  swing  round  towards  the 

point  of  minimum  density,  as  predicted  by  the  theory.  After  July  4,  however,  the  drag 

takes  command  and  $  rapidly  approaches  <p  ,  nearly  coinciding  with  it  at  the  end 
P  P* 

(the  difference  being  within  the  errors  caused  by  assuming  an  equatorial  position  for  the 
maximum  and  minimum  densities). 


The  values  of  zf  and  up  have  also  been  calculated  for  each  day  of  the  orbit 
determination,  and  they  are  plotted  as  triangles  in  the  (z  cos  to,  z  sin  w)  plane  in 
Fig  13.  The  corresponding  observational  values  of  z  cos  to  and  z  sin  to  are  plotted 
as  circles,  the  numbers  indicating  the  orbit  number  in  Tables  1  and  2. 


Fig  13  can  be  interpreted  satisfactorily  in  the  light  of  the  theory,  which  suggests 
that  che  observational  values  should  be  moving  towards  the  appropriate  triangles  as  soon 


enough  to  overcome  che  effects  of  the  gravitational  field.  For 
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orbits  1  to  5  the  initial  direction  of  motion  under  the  influence  of  the  gravitational 

field  is  still  maintained,  though  it  also  happens  to  be  in  the  direction  required  by  the 

minimum-density  criterion,  ie  the  values  shown  as  circles  are  heading  towards  the 

appropriate  triangles.  For  orbits  6  to  10,  when  the  day-to-night  variation  in  density 

exerts  more  influence,  the  observed  values  of  z  head  straight  towards  the  appropriate 

Zp  values.  For  orbits  10  to  12,  z?  moves  rapidly  to  zero  and  out  again,  and  the 

observed  values  of  z  head  towards  the  values  of  z?  ,  ignoring  the  short-lived 

excursion  of  z  to  zero.  After  orbit  13  the  height  of  the  satellite  has  decreased  to 

r 

below  190  km  where  the  day-to-night  variation  is  much  smaller,  with  F  <  0.1  ,  so  the 
tendency  for  z  to  approach  z?  grows  slightly  weaker,  while  there  is  an  increasing 

tendency  for  e  to  decrease  due  to  the  basic  spherical-atmosphere  drag. 

At  this  stage  it  becomes  preferable  to  interpret  the  results  in  terms  of  e  rather 
than  z  .  So  far  i£  has  been  possible  to  treat  z  (»  ae/H)  and  e  as  similar,  because 
the  change  in  H  between  one  orbit  and  the  next  has  not  been  large  enough  to  invalidate 
the  assumption  of  constant  H  in  the  theory.  But  for  orbits  12  to  15  there  are 
substantial  changes  in  H  ,  and  so  it  is  better  to  consider  ef  -  Hz^/a  (where  H/a  has 
separate  values  for  each  orbit)  and  to  look  at  the  variation  of  e  relative  to  ep  . 

The  values  of  e  and  e?  on  orbits  12  to  15  are  shown  on  the  inset  diagram  in  Fig  13. 

In  the  last  three  days  the  value  of  slowly  increases  from  335°  to  350°.  The 

actual  value  of  u  is  rapidly  decreasing  between  orbit  12  and  orbit  13,  but  then  slows 
down  between  orbits  13  and  14  (see  Fig  11)  as  it  has  overshot  the  value  of  i»F  ,  and 
increases  between  orbits  14  and  15.  On  orbits  12  and  13,  e  <  e?  and  so  e  would  be 
expected  to  increase,  and  does  so.  Between  orbits  14  and  15,  ep  becomes  less  than  e  , 
and  the  slight  decrease  in  e  on  orbit  15  reflects  this.  Thus  the  theory  is  fully  borne 
out  by  the  observed  variation. 

5.3  Analysis  of  perigee  position 

Having  established  that  the  movement  of  perigee  is  in  conformity  with  the  theory, 
it  is  worth  looking  at  the  motion  of  the  perigee  point  in  more  detail,  to  examine  its 
progress  towards  the  minimum-density  point,  while  remaining  within  the  orbital  plane. 

Fig  14  gives  sketches  showing  how  this  progress  develops.  Fig  14a  defines  the 
notation:  P  is  the  actual  perigee  position;  D  is  the  point  of  minimum  density 
(assuming  an  equatorial  Sun,  as  in  the  theory);  and  Pj,  is  the  'final'  value  of  perigee 
towards  which  perigee  should  be  moving  according  to  the  theory,  te  the  point  of  minimum 
density  in  the  orbital  plane. 

Fig  14b-e  illustrate  the  movement  of  P  on  four  selected  orbits.  Nos.  1,4,  10 
and  14.  Initially,  Fig  14b,  P  is  near  northern  apex  (w  *  120°)  and  <u  is  slowly 
increasing  due  to  the  gravitational  field.  At  this  time  P  is  at  local  time  16  h.  By 
orbit  4  the  argument  of  perigee  has  ceased  to  increase  and  the  effect  of  the  day-to-night 
variation  in  air  density  gradually  begins  to  take  command  over  the  effect  of  the 
gravitational  field.  However,  the  perigee  happens  to  be  close  to  the  maximum-density 
point  on  orbits  1  to  4,  and  so  at  this  stage  P  has  a  long  way  to  go  to  reach  Pp  .  In 
orbit  10,  shown  in  Fig  14d,  P  is  rapidly  moving  towards  Pf  ,  the  minimum-density  point 
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in  the  orbital  plane.  Finally,  in  Fig  I4e,  a  sketch  of  orbit  14,  the  perigee  position 
has  slightly  overshot  the  position  of  Pp  .  (However,  as  stated  previously,  the  angular 
separation  of  P  and  Pp  is  within  the  errors  caused  by  assuming  an  equatorial  position 
for  the  Sun.) 

6  THE  REMAINING  ORBITAL  ELEMENTS 

6. 1  The  semi  major  axis 

The  values  of  semi  major  axis  a  from  Tables  1  and  2  are  plotted  as  circles  in 
Fig  IS  and  joined  by  a  smooth  curve  (full  line).  They  are  plotted  against  t'/tL  where 
t'  is  the  time  in  days  after  June  28.0  and  t^  is  the  lifetime  from  June  28.0, 

■ie  13.69  days. 

The  dashed  curve  in  Fig  15,  determined  from  equation  (4.90)  of  Ref  18,  gives  the 
theoretical  change  in  a  over  the  same  interval  as  the  orbit  determinations,  assuming 
no  variation  of  density  with  time  and  a  constant  density  scale  height  H  .  If  a^  is 
the  initial  value  of  a  (6624.9  km),  equation  (4.90)  may  be  written 


(ID 


where  n  ■  1  -  exp|-  (aQ  -  a^)/H|  and  \  ,  taken  as  6472  km,  is  the  value  of  a  at 

decay.  Equation  (II)  applies  for  a  spherical  or  oblate  atmosphere. 

The  dash-dot  curve  gives  a  nearer  approximation  by  allowing  for  a  variable  scale 
height,  as  in  equation  (6.81)  of  Ref  18,  where  the  symbol  u  represents  the  rate  of 
increase  of  H  with  height  and  the  variation  of  a  with  t'  is  given  by 


(a0  ~  a)  ( 
H8  i 


iuU  + 


0(iii2) 
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Here  Hg  is  the  value  of  H  at  one  scale  height  below  the  initial  height  ( ie  at  a 
distance  ag  -  Hg  from  the  Earth’s  centre),  n  is  calculated  with  H  ■  Hg  , 

(n  *  0.986  here),  and  U  is  a  function  of  t'/tL  plotted  in  Fig  6.5  of  Ref  18. 

In  the  theory  it  is  assumed  that  y  <  0.2  .  Here  CIRA  1972  indicates  that  y  is 

0.2  on  June  28  and  increases  to  0.35  on  July  11,  and  a  2Z  error  in  t'  is  introduced  by 

,  2 
limiting  y  to  0.2.  There  is  also  an  error  due  to  the  neglect  of  terms  of  order  Jy 

in  equation  (12),  which  is  0.02  if  y  -  0.2. 

Even  with  the  limitations  of  equation  (12),  however,  the  dash-dot  curve  gives  a 

better  approximation  to  the  true  decrease  in  semi  major  axis  than  the  constant-H  equation 

(dash  curve).  The  maximum  departure  of  the  actual  t'/t^  from  the  theoretical  value  Is 

4Z  with  equation  (11)  and  2 X  with  equation  (12).  This  is  well  within  the  expected  error 
2 

due  to  neglect  of  y  and  variations  of  density  with  time. 

6.2  The  right  ascension  of  the  ascending  node 

The  values  of  right  ascension  of  the  ascending  node,  ft  ,  are  plotted  in  Fig  16 
against  date  and  they  follow  the  expected  pattern,  decreasing  by  about  5.7  deg/day  due  to 
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the  effect  of  the  even  zonal  harmonics  in  the  geopotential.  If  there  was  any  variation  in 
fi  due  to  atmospheric  rotation,  the  perturbations  due  to  the  even  zonal  harmonics  and 
lunisolar  forces  would  first  have  to  be  removed,  and  then  the  residual  variation  fitted 


(13) 


Here,  I2/Iq  <  because  e  <  0.001  and  z  <  0.136  ,  and  the  effect  would  be 

negligible.  So  it  is  not  possible  to  determine  A  from  the  variation  of  0  . 

6.3  The  orbital  decay  rate 

Fig  17  shows  as  circles  the  orbital  decay  rates  n  (■  2M2)  in  deg/day2  for  the 
orbits  of  Tables  1  and  2,  the  values  for  the  orbits  13A  and  B  and  14A  and  B  being  plotted 
at  times  away  from  epoch,  for  the  reasons  given  in  section  2.  These  values  of  A  could 
have  been  used  to  determine  values  of  density,  but  they  are  instantaneous  values,  and  the 
method  preferred  (section  4)  was  to  use  the  average  daily  change  in  semi  major  axis  Aa 
to  give  daily  values  of  a  (*  Aa/At)  .  These  values  of  Aa  are  also  plotted  in  histogram 
form  in  Fig  17. 

Although  the  two  sets  are  not  strictly  comparable,  they  should  of  course  show  the 

same  trends.  The  means  of  the  successive  daily  values  of  Aa  ,  which  represent  the  mean 

drag  over  2  days,  are  joined  with  a  broken  line  up  to  July  9  and  this  closely  parallels 

2 

the  unbroken  curve.  The  ratio  of  2M2(deg/day  )/Aa(km)  should  be  1.31  on  28  June, 
increasing  to  1.34  on  9  July;  but  after  July  9,  with  decay  imminent,  the  ratio  of  an 
instantaneous  to  an  average  value  obviously  becomes  an  unacceptable  comparison. 

7  DISCUSSION  AND  CONCLUSIONS 

7. 1  The  orbit  determinations 

The  orbit  has  been  determined  daily  from  June  28  to  July  1 1  and  the  sets  of 
elements  are  given  in  Tables  1  and  2.  In  Table  2  three  sets  of  elements  for  epochs 
July  10.0  and  11.0  are  given,  but  the  orbit  C  for  both  epochs  is  reconmended  because  it 
was  determined  from  observations  spanning  the  epoch.  The  elements  for  orbit  15  (July 
11.5)  in  Table  2  are  converted  from  the  orbital  elements  (orbit  14D)  derived  from  the 
fitting  of  the  observations  between  08.08  and  16.02  UT  on  July  11. 

The  standard  deviations  of  the  orbital  elements  in  the  daily  orbits  from  June  28.0 
to  July  11.0  correspond  to  radial  and  cross-track  accuracies  of  30  to  40  m.  If  these 
accuracies  are  realistic,  the  orbits  are  the  most  accurate  ever  published  for  such  a 
high-drag  satellite.  Such  a  conclusion  would  not  be  too  surprising,  because  the  satellite 
was  intensively  observed  by  all  the  radars  of  the  North  American  Air  Defense  Cotmand,  and 
also  by  the  US  Navy's  Navspasur  system  and  by  several  French  radars,  and  all  these 
observations  were  successfully  used  in  the  orbit  determination.  Furthermore  the  orbit 
refinement  program  used,  Gooding's  PROP  6,  has  proved  itself  to  be  extremely  reliable  for 
high-drag  orbits,  whenever  the  observations  are  confined  to  a  time  span  of  a  day  or  two 
and  the  irregular  variations  in  air  drag  are  not  too  severe.  Both  these  criteria  are 
satisfied  here,  and  there  is  no  reason  to  suppose  that  errors  of  more  than  20  m  radial  or 
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cross-track  arise  through  errors  in  the  PROP  model.  Despite  these  advantages,  an 
accuracy  of  30  to  40  m  does  seem  somewhat  over-optimistic,  (a)  because  the  observations 
have  a  basic  accuracy  of  100  to  200  m,  and  (b)  because  Skylab  had  a  span  of  27  m  and  a 
length  of  26  m,  and  reflections  could  have  come  from  different  parts  of  the  spacecraft. 

All  in  all,  an  accuracy  of  about  60  m  radial  and  cross-track  seems  more  realistic,  and 
this  is  borne  out  by  the  analysis  of  inclination  (section  3). 

The  fitting  of  the  observations  between  08.08  and  16.02  UT  on  July  1 1  -  up  to 
within  35  minutes  of  the  burn-up  -  was  remarkably  good,  and  there  seems  no  reason  why  the 
accuracy  of  this  orbit  (radial  or  cross-track)  should  be  very  much  worse  than  the  others. 
Degradation  by  a  factor  of  about  2  might  be  expected  and  this  is  in  conformity  with  the 
standard  deviations  obtained.  So  the  orbit  for  noon  on  the  last  day  should  be  accurate  to 
about  100  m  radial  and  cross-track,  over  the  8-hour  time  interval  of  the  observations. 

7.2  Other  orbital  work  on  Skylab  1 

The  decay  of  Skylab  generated  great  interest  because  it  was  known  that  large  pieces 
would  survive  decay,  and  these  were  seen  as  a  threat  to  densely  populated  areas.  NORAD/ 

NASA  were  issuing  updated  orbits  at  frequent  intervals  and  these  have  been  collected  and 
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presented  by  Wakker  from  April  1978  until  re-entry  in  July  1979.  Carrou  records  how 

the  French  covered  the  decay  and  gives  plots  of  semi  major  axis  and  inclination  from 

April  to  July  1979.  However  the  main  purpose  of  these  determinations  was  to  predict 

re-entry  and  not  to  analyse  the  orbital  elements,  and  the  cross-track  accuracy  achieved 

appears  to  be  about  300  m. 

.  .  27 

The  orbit  of  Skylab  1  was  determined  by  Brookes  and  Moore  at  46  epochs  between 

January  1974  and  August  1976  using  PROP  6  with  2100  observations  from  five  sources;  the 

Malvern  Hewitt  camera,  the  Cape  kine theodolite,  the  US  Navy  sensors,  the  Finnish 

theodolite  and  visual  observers.  Their  orbits,  which  were  used  in  determining  air 

density,  had  average  accuracies  of  75  m  and  130  m  radial  and  cross-track  respectively, 

,  2 
but  were  of  course  at  a  time  when  drag  was  slight  (M_  averaged  0.01  deg/day  ,  as  compared 

2  * 
with  1.6  to  60  deg/day  here). 

7.3  The  analysis  of  the  orbits 

The  values  of  inclination  have  been  analysed  to  obtain  an  average  atmospheric 

rotation  rate,  A  .  The  value  found  (Fig  6)  was  1.10  ±  0.07  for  an  average  height  band  of 

210  to  220  km.  This  result,  obtained  after  removing  the  resonance  perturbations,  conforms 

well  with  the  results  in  Ref  12,  where  a  value  of  1.1  for  A  is  indicated  for  average 

conditions.  The  values  of  inclination  at  16th-order  resonance  have  been  successfully 

analysed  (Fig  5),  and  lumped  16th-order  geopotential  coefficients  C*  and  S.!  have 

lb  lb 

been  derived  for  the  first  time  from  resonance  analysis.  The  values  obtained  were 
9-0,1  o_0,l 

I0’c,'  -  147  ±  42  and  ICTS^  -  4  ±  35  . 

The  air  density  has  been  determined  at  daily  intervals  at  appropriate  daily 
values  of  height.  The  density  values  have  been  standardized  to  a  fixed  value  of 
exospheric  temperature  and  are  compared  in  Fig  8  with  values  of  density  derived  from 
CIRA  1972.  For  the  first  8  days  the  results  indicate  densities  I5Z  lower  than  the 
annual-averaged  CIRA  values;  however,  the  CIRA  model  gives  a  decrease  of  14Z  in  early 


July  due  to  the  semi-annual  variation.  So  it  appears  that  the  semi-annual  variation  in 

July  1979  conformed  to  that  specified  by  the  CIRA  1972  model,  in  contrast  to  results  from 

2 1 

the  early  1970s,  when  much  larger  semi-annual  variations  occurred  . 

The  variation  in  eccentricity  for  Skylab  has  been  analysed  in  detail  and  is  found 
to  follow  that  specified  by  the  theory  for  an  atmosphere  with  day-to-night  variation  in 
air  density.  In  particular  the  observed  increase  in  eccentricity  a  few  days  before  decay 
(Fig  9)  can  be  fully  explained  by  the  theory. 

It  has  been  shown  that,  in  accordance  with  the  theory  for  near-circular  orbits  about 
to  decay  in  an  atmosphere  with  day-to-night  variation  in  density,  the  perigee  moves  towards 
the  point  of  minimum  density  in  the  orbit.  This  process  is  illustrated  in  Fig  14.  The  com 
bined  movement  of  eccentricity  and  argument  of  perigee  also  conforms  to  theoretical  expecta 
tions,  as  shown  in  Fig  13. 
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Appendix 

MEAN  HEIGHT  AT  WHICH  VALUES  OF  p  OR  A  APPLY 


WHEN  DETERMINED  FROM  CIRCULAR  ORBIT  DECAY 


If  air  density  is  determined  from  the  decay  rate  of  a  circular  orbit  with  radius 
that  decreases  from  a^  initially  to  a^  ,  at  what  value  a  of  a  does  the  density 
apply?  To  a  first  approximation,  1  -  j(ag  +  aj)  ,  but  a  more  accurate  value  is  required 
here. 

From  equation  (4.84)  of  Ref  18,  the  change  A* a  in  one  revolution  due  to  drag  in  a 
spherically  symmetrical  atmosphere  is 

A' a  *  -  2ir6a^p  ,  (A—  1 ) 


where  6  is  the  area/mass  parameter  and  p  is  the  air  density  at  distance  a  from  the 
Earth's  centre,  given  by  .  . 


p  -  pQ  exp|g(a0  -  a) |  . 


In  equation  (A-2),  p^  is  the  density  at  distance  aQ  and  1 / B  “  H  is  the  density 
scale  height.  The  rate  of  decrease  of  a  at  any  time,  da/dt  ,  is  A'a/T  ,  where  T  is 
the  orbital  period,  equal  to  2ir/n  ,  where  n  is  the  mean  motion.  Thus,  from  (A-l), 


da  .2 

*  -  n6a  p  , 


(A-3) 


whence 


exp{S(a  -  aQ)^f 


-  n6a  p. 


from  (A-2).  Since  na  does  not  vary  by  more  than  0.2%  in  one  day  for  1973-27A,  we  may 
take 

nSa^Pg  -  k  (A-5) 

as  constant,  so  that  (A-4)  may  be  integrated  between  t  *  tg  and  t  -  t(  to  give 


1  -  exp^SUj  -  aQ)|  -  8k(tj  -  tQ)  . 


Since  the  rate  of  change  of  a  is  proportional  to  p  by  (A-3),  the  mean  density 
p  befveen  times  t^  and  tj  is  given  by 


_x_  /Bi£  . 


on  using  (A-3).  From  (A-7)  and  (A-5). 
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Appendix 


On  using  (A-2)  and  (A-6),  equation  (A-8)  becomes 


which  gives 


(A-9) 


If  6(a0 
give 


aj)  is  small,  equation  (A-9)  can  be  expanded  in  powers  of  B(aQ  -  a^  to 


(A-10) 


For  1973-27A,  equation  (A-10)  is  adequate  except  for  the  last  day,  when  B (a^  -  a^)  *  0.82 
and  (A-9)  is  needed. 

2 

The  effects  of  atmospheric  oblateness  introduce  a  factor  (1  +  ic  )  into  the  right- 
hand  side  of  equation  (A-l),  as  shown  by  equation  (5.49)  of  Ref  18.  Since  this  differs 
from  1  only  by  second-order  terms,  and  is  constant,  the  analysis  should  also  be  valid  for 
an  oblate  atmosphere. 

Since  the  change  in  inclination  Ai  is  proportional  to  p  ,  the  same  analysis 
should  be  valid,  and  the  values  of  rotation  rate  A  obtained  by  fitting  the  decrease  in 
inclination  may  be  taken  as  applying  at  the  distance  a  given  by  (A-9). 
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Fig  4  Values  of  inclination  from  Fig  3  plotted  against  period,  with  fitted  straight  line 
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resonance  with  the  geopotential  have  been  successfully  analysed,  to  give  the  first 
values  of  16th-order  geopotential  coefficients  determined  from  resonance,  and  a 
value  of  1.10  ±  0.07  rev/day  for  the  atmospheric  rotation  rate  at  a  height  of 
210-220  km. 

The  daily  changes  in  semi -major  axis  have  been  used  to  determine  13  daily 
values  of  air  density,  at  heights  from  252  km  down  to  179  km.  All  agree  well  with 
the  CIRA  1972  model,  and  indicate  a  smaller  semi-annual  variation  than  in  the  early 
1970s. 

The  variations  of  eccentricity  and  argument  of  perigee  take  unusual  forms,  but 
detailed  analysis  shows  that  the  variations  are  in  full  accord  with  the  theory  for 
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